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Abstract.  We  consider  a  model  for  the  interrogation  of  a  planar  Debye  medium  by  a 
non-harmonic  microwave  pulse  from  an  antenna  source  in  free  space,  and  we  compute  the 
reflected  solution  using  finite  elements  in  the  spatial  variables  and  finite  differences  in  the 
time  variable.  Perfectly  Matched  Layers  (PMLs)  and  an  absorbing  boundary  condition  are 
used  to  damp  waves  interacting  with  artificial  boundaries  imposed  to  allow  finite  compu¬ 
tational  domains.  We  present  simulation  results  showing  that  numerical  reflections  from 
interfaces  at  PML  boundaries  can  be  controlled. 


1.  Introduction 

The  purpose  of  this  paper  is  to  demonstrate  computationally  that  one  can  implement 
a  two-dimensional  version  of  the  electromagnetic  interrogation  problems  introduced  in  [3]. 
Here  we  use  perfectly  matched  layers  (PMLs)  as  absorbing  layers  at  artificial  boundaries  used 
to  define  finite  computational  domains.  We  carry  out  calculations  to  verify  that  artificial 
reflections  do  not  contaminate  reflections  from  dielectric  layer  interfaces  used  to  determine 
dielectric  parameters  as  well  as  physical  geometry  in  inverse  problem  formulations.  Our 
results  provide  evidence  that  the  one- dimensional  normally  incident  plane  wave  ideas  in  [3] 
can  be  extended  to  treat  higher  dimensional  problems  in  which  obliquely  incident  waves  play 
a  nontrivial  role. 

In  [3]  the  authors  develop  a  theoretical  and  computational  framework  to  use  pulsed  electro¬ 
magnetic  interrogating  signals  to  determine  dielectric  and  geometric  properties  of  materials. 
This  framework  involves  time  domain  computations  of  electromagnetic  signals  from  an  an¬ 
tenna  through  vacuum  to  the  target  and  return.  Even  in  the  one-dimensional  case  treated 
in  [3]  (where  one  uses  an  infinite  antenna  sheet  and  polarized  plane  waves  to  achieve  a  one¬ 
dimensional  finite  spatial  domain),  the  problems  are  computationally  intensive.  Moreover, 
there  are  several  difficulties  in  extending  this  methodology  to  two  and  three  dimensions  in 
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addition  to  the  usual  increased  computational  complexities  involved  in  moving  to  higher 
spatial  dimensions.  First,  interrogating  signals  from  a  finite  antenna  will  produce  oblique 
incident  waves  to  the  target  and  these  must  be  treated  in  the  reflections.  The  uniformity 
assumptions  made  in  [3]  to  yield  one-dimensional  finite  spatial  domains  will  not  be  applica¬ 
ble  and  an  infinite  spatial  domain  must  be  approximated  by  a  finite  computational  domain 
with  artificial  boundaries.  Since  perfectly  absorbing  boundary  conditions  are  not  available 
in  higher  dimensions,  some  type  of  boundary  damping  must  be  formulated  so  that  artificial 
reflections  will  not  interfere  with  reflections  from  the  target. 

In  the  treatment  here,  we  model  the  propagation  of  a  non-harmonic  pulse  from  a  finite 
antenna  source  in  free  space  across  a  planar  interface  into  a  dielectric.  The  dielectric  is  a 
Debye  medium  with  Ohmic  conductivity.  We  use  finite  elements  in  the  spatial  variables 
and  finite  differences  in  the  time  variable  to  compute  the  components  of  the  electric  held  in 
the  case  where  the  signal  and  dielectric  parameters  are  independent  of  the  y  variable  (the 
only  uniformity  assumption  made  here).  Figure  1  depicts  the  antenna  and  dielectric  slab 
geometry  we  use  in  our  problem  formulation  with  the  infinite  dielectric  slab  perpendicular 
to  the  2-axis  and  uniform  in  the  region  z\  <  z  <  Z2-  The  uniform  strip  antenna  is  located 
in  the  xy- plane,  is  infinite  in  the  ^/-direction  and  finite  in  the  x- direct  ion  lying  in  the  region 
— oo  <  y  <  oo,  —x  <  x  <  x.  An  alternating  current  along  the  x- direction  then  produces 
an  electric  held  that  is  uniform  in  y  with  nontrivial  components  Ex  and  Ez  depending  on 
( t ,  x,  z)  which,  when  propagated  in  the  xx-plane,  results  in  oblique  incident  waves  on  the 
dielectric  surface  in  the  xy-plane  at  z  =  Z\. 

We  also  assume  that  the  dielectric  is  backed  by  a  supraconductive  material  with  an  infi¬ 
nite  conductivity,  so  one  side  of  the  computational  domain  will  have  a  perfectly  reflecting 
boundary  condition.  Artificial  boundaries  on  the  other  three  sides  are  assumed,  producing 
an  approximating  finite  computational  domain.  Energy  will  also  rehect  oh  the  boundaries  on 
these  other  three  sides  of  the  computational  domain,  and  there  is  a  critical  need  to  prevent 
or  delay  this  energy’s  return  to  the  domain  of  interest.  We  use  Perfectly  Matched  Layers 
to  absorb  incident  waves  without  reflection  and  damp  the  absorbed  waves,  and  we  use  a 
(partially)  absorbing  boundary  condition  to  limit  the  amount  of  energy  that  reflects  oh  a 
boundary  of  the  computational  domain.  When  the  outgoing  wave  returns  to  the  domain 
of  interest  after  twice  traversing  the  PMLs,  the  wave  can  be  sufficiently  attenuated  to  be 
negligible  when  compared  to  the  rehections  from  the  dielectric. 

Berenger  [5]  formulated  construction  of  PMLs  for  free  space  using  held  splitting.  Sub¬ 
sequently  Sacks  et.  al.  [13]  developed  the  anisotropic  formulation  for  PMLs  that  we  follow 
here.  We  formulate  the  PMLs  in  terms  of  auxiliary  equations  similar  to  the  second  order 
Lorentz  equations  employed  in  [3].  (A  similar  auxiliary  equation  PML  formulation  can  also 
be  found  in  [14]  and  [15].)  The  absorbing  boundary  condition  we  use  for  PMLs  is  adapted 
from  Engquist  and  Majda’s  hrst  order  absorbing  boundary  condition  for  free  space  [10]. 

A  PML  can  be  constructed  to  match  impedance  with  either  free  space  or  a  dielectric. 
Each  coordinate  axis  will  determine  a  one-parameter  family  of  PMLs.  The  PML  parameter 
determines  the  degree  of  energy  attenuation  in  the  PML.  The  parameter  domain  is  [0,  oo) , 
and  the  PML  with  parameter  0  corresponds  to  free  space  or  the  dielectric  depending  on 
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Figure  1.  Antenna  and  dielectric  material  geometry 

whether  the  impedance  of  the  one-parameter  PML  family  is  matched  to  free  space  or  the 
dielectric.  All  planar  interfaces  between  two  PMLs  from  the  same  family  will  be  reflectionless 
if  the  interface  is  normal  to  the  family’s  coordinate  axis.  Using  one  or  multiple  PMLs  from  the 
appropriate  family,  one  can  construct  a  buffer  region  to  a  vacuum  or  to  a  dielectric  that  will 
absorb  energy  without  reflection  and  then  attenuate  that  energy.  In  numerical  computation 
there  will  be  reflections  at  the  interfaces  between  PMLs  from  the  same  family,  but  these 
reflections  can  be  controlled  by  making  the  jumps  in  the  PML  parameter  sufficiently  small. 
We  note  that  there  are  relatively  few  efforts  in  the  research  literature  on  using  PMLs  with 
finite  elements  (see  [7]  for  recent  efforts  and  related  references)  and,  more  specifically,  with 
finite  antenna  generated  time  domain  Maxwell  equation  propagated  fields. 

Our  research  extends  the  work  of  Banks,  Buksas,  and  Lin  [3]  who  use  finite  elements 
in  space  and  finite  differences  in  time  to  solve  a  corresponding  one-dimensional  forward 
scattering  problem.  However,  in  one  dimension,  PMLs  are  not  needed  since  an  absorbing 
boundary  condition  will  absorb  all  incident  energy.  For  example,  if  waves  propagate  at  unit 
speed,  the  boundary  condition  £  —  ^  at  z  —  Zo  absorbs  all  waves  moving  in  the  negative 
^-direction.  There  is  no  reflection  since  a  reflected  wave  would  have  the  form  f(z  —  t ). 
Unfortunately,  absorbing  boundary  conditions  in  higher  space  dimensions  will  reflect  some 
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energy  from  all  plane  waves  at  non-normal  incidence  [10].  In  our  formulation,  PMLs  provide 
the  necessary  extra  damping  to  attenuate  the  energy  that  escapes  the  absorbing  boundary 
condition. 

As  we  have  noted  above,  extending  the  work  in  [3]  to  two  dimensions  also  changes  the 
incidence  of  the  electric  field  on  the  vacuum-dielectric  interface  from  normal  to  oblique 
incidence.  Blaschak  and  Franzen  [6]  use  Fourier  series  in  the  frequency  domain  to  compute 
the  propagation  of  a  time  harmonic  pulse  train  of  plane  waves  that  enter  a  dielectric  across  a 
planar  boundary  at  oblique  incidence.  They  showed  that  precursors  are  excited  by  short-rise- 
time  pulses  at  oblique  incidence.  However,  the  use  of  Fourier  series  restricts  this  approach 
to  harmonic  pulses. 

We  begin  with  Maxwell’s  equations,  introduce  constitutive  relations  for  a  diagonally 
anisotropic  material,  and  derive  the  wave  equation  and  its  corresponding  variational  formula¬ 
tion  (Section  2).  We  then  compute  the  reflection  coefficient  at  a  planar  interface  between  two 
diagonally  anisotropic  materials  (Section  3)  and  find  sufficient  conditions  for  a  reflectionless 
interface. 

In  Section  4  we  formulate  our  problem  geometry  and  the  parametrized  families  of  PMLs 
which  we  use.  We  give  a  single  variational  equation  for  the  electric  field  that  is  valid  on  the 
entire  computational  domain,  and  we  identity  differential  equations  that  are  satisfied  by  the 
polarization  and  PML  terms  in  the  variational  equation.  A  perfectly  reflecting  boundary 
condition  is  applied  to  the  back  of  the  dielectric,  and  an  absorbing  boundary  condition  is 
applied  to  the  other  three  sides  of  the  computational  domain. 

We  then  discretize  our  computational  domain  and  derive  a  finite  dimensional  system  of 
equations  for  the  electric  field  (Section  5).  We  give  numerical  results  in  Section  6  which  verify 
that  PMLs  significantly  attenuate  energy  and  that  numerical  reflections  from  interfaces  at 
PML  boundaries  can  be  controlled. 

2.  The  Wave  Equation 

To  derive  the  wave  equation  for  a  diagonally  anisotropic  region,  we  begin  with  Maxwell’s 
equations  for  a  region  with  free  charge  density  p  =  0: 

V  ■  D  =  0 

V  ■  B  =  0 

(1) 

V  x  E  =  -dtB 

V  x  H  =  dtD  +  J. 

The  vectors  in  (1)  are  functions  of  position  r  =  (x,y,z)  and  time  t.  The  current  density 
is  the  sum  of  a  conduction  current  density  and  a  source  current  density:  J  =  Jc  +  Js.  We 
assume  only  free  space  can  have  a  source  current  (e.g.,  in  our  antenna),  and  thus  either 
Jc  =  0  or  Js  =  0,  depending  on  whether  or  not  the  region  is  free  space. 

We  next  introduce  constitutive  relations  that  are  sufficiently  general  to  allow  for  diagonal 
anisotropy.  Let  iS+  be  the  space  of  tempered  distributions  on  (— oo,  oo)  supported  in  [0,  oo). 
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With  convolution  as  multiplication,  5+  is  a  commutative  ring  with  identity  ([11],  Sections 
5.3,  8.3),  and  the  Dirac  delta  function  6  is  the  multiplicative  identity.  Let  _Ad+  be  the 
multiplicative  group  of  invertible  elements  in  <S+,  let  C[0,oo)  be  the  space  of  continuous 
functions  on  [0,  oo),  and  define 

(2)  B+  =  {L  :  L  is  a  Injection  on  C[ 0,  oo)  and  Lf  =  g  *  f  for  some  g  G  M.+}. 

Observe  that  B+  is  a  group  when  the  binary  operator  is  composition  of  mappings.  In 
particular  L  G  B+  if  and  only  if  L_1  G  B+ . 

Let 

T  =  diag  (Tx,Ty,Tz) 

where  Tx,  T  ,  Tz  G  B+,  and  let  the  diagonally  anisotropic  material’s  constitutive  relations 
be  given  by 

D  =  g  * ( TE ) 

(3)  B  =  g0TH 

Jc  =  aTE 

where  g  G  .Ad+,  convolution  is  in  the  time  variable,  /i()  is  the  permeability  of  free  space,  and 
the  constant  a  >  0  is  the  conductivity.  These  constitutive  relations  are  sufficiently  general  to 
include  free  space  or  an  isotropic  dielectric  with  Debye  model  (as  well  as  many  other  models) 
for  polarization.  For  instance,  if  T  is  the  identity  operator,  a  =  0,  and  g  =  eo<5  where  eo  is 
the  permittivity  of  free  space,  then  (3)  describes  a  vacuum. 

Convolution  in  time  commutes  with  the  divergence  operator: 

V  •  D  =  V  •  (g  *  TE) 

=  g*{V-TE). 

Since  g  G  M.+  and  V  •  D  —  0  it  follows  that  V  ■  TE  =  0.  Let  TL(t)  be  the  Heaviside  function. 
Using  the  constitutive  relations  (3),  we  can  rewrite  Maxwell’s  equations  (1)  as 

V  ■  TE  =  0 

V  •  TH  =  0 

(4) 

V  x  E  =  -fioT(dtH) 

VxH  =  T(dt((g  +  gH)  *  E))  +  Js. 


t  =  pyrpw'r 

ft  =  (T„TyT,)-% 


t  =  x,y,z. 


Then  T  =  diag(T1:,  T  Tz 


6 


H.  T.  BANKS  AND  BRIAN  L.  BROWNING 


We  next  derive  the  wave  equation  for  a  region  where  T  and  g  are  constant  in  the  spatial 
variables.  Since  T  will  also  be  spatially  constant  in  Q,  V  •  TE  —  0  by  (4).  A  calculation 
shows  that 

-(V  .  (TV))E  =  -(V  .  (TV))E  +  V(V  ■  (TE)) 

(o) 

=  T-\V  x  (T-1(V  x  E))). 

Consequently  using  (6)  and  (4)  we  obtain 

(-V  •  (TV))E  =  T-1(V  x  T_1(V  x  E)) 

=  -g0dtT~\\7  x  H) 

=  -g0  d2t({g  +  oH)  *  E)  -  goT~ldJs 
=  -fi0dt(g  *  E)  -  Hq(j dtE  -  HQ T~ldtJs. 

Thus  the  wave  equation  for  E  in  a  homogeneous,  diagonally  anisotropic  region  with  con¬ 
stitutive  relations  (3)  is  given  by 

(7)  no  d2t{g  *  E)  +  ^adtE  -  (V  •  (TW))E  =  -/r0  T~ldtJs. 


Dehne  the  electric  polarization  vector  P  by 

g  *  E  =  Co 6rE  T  P 

where  er  is  a  relative  permittivity  that  incorporates  the  instantaneous  polarization  [3,  p.  11] . 
Then  the  wave  equation  (7)  becomes 

(8)  /doeoerE  +  go P  +  go^E  —  (V  •  {TW))E  =  —goT  1  Js 

where  '  =  dt  and  ”=  <9t2.  In  our  problem  a  current  source  will  exist  only  in  the  antenna  where 

T  is  the  identity  operator,  so  the  right  side  of  (7)  and  (8)  will  reduce  to  either  g0Js  or  0. 

The  components  of  (8)  can  be  expressed  in  variational  form.  Let  have  piecewise  smooth 
boundary,  and  let  0  be  a  test  function  on  Q.  For  a  function  /  on  dehne 

</»=  [  mm  dv. 

Jn 

If  /  and  (j)  are  also  defined  on  dtt  let 

(/,  4>)du  =  [  f{r)<t>(f)  ds. 

J  dfl 

We  are  interested  in  computing  the  ^-component  of  E  which  is  the  signal  measured  by 
the  finite  antenna  (see  Figure  1)  when  reflected  by  the  dielectric.  Let  x  denote  the  standard 
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basis  vector  parallel  to  the  x-axis,  and  define 

(9) 


E  =  E  ■  x 
P  =  P  ■  x 
Js  =  (T”1  Ja)  ■  x. 


Applying  the  product  rule,  the  divergence  theorem  and  the  definition  of  T  in  (5),  we  obtain 

((V  •  {TV))E,  (/>)  =  [  V  •  {<t>T VE)  dV  -  f  (V0)  •  (TVE)  dV 

Jn  Jn 

=  [  4>(T'S7E)  •  rids  -  [  (V</>)  •  (TVE)  dV 
Jon  Jn 

=  <$(? [Elmfion-  (diiTMdit) 


i=x,y,z 


i=x,y,z 


where  n  =  (nx,ny,nz)  is  the  outward  unit  normal  to  the  boundary  <9fh  Consequently  the 
variational  form  of  the  wave  equation  (8)  on  is 

(10)  n060(erE,(/))  +/x0(P»  +n0{crE,(f))  +  ^  (dpT^E),  8$)  -  ^  {di(TiE),ni<j))dn 


i=x,y,z 


i=x,y,z 


ho  (JsiP)  i 


or  more  succinctly 

H0e0(6rE,<f))  +  /xo(P,0)  +  /J,0(aE,  (/>)  +  (VTA,  V</>)  -  (VT^,n0)aa  =  -/x0(  js,  </>). 

Our  computational  domain  will  be  a  region  T>  that  can  be  partitioned  into  subregions 
{Qfc}  with  piecewise  smooth  boundaries  such  that  T  and  g  are  constant  in  the  spatial  vari¬ 
ables  on  each  O*,.  Consequently,  the  components  of  a  solution  of  Maxwell’s  equations  (1) 
with  constitutive  relations  (3)  will  satisfy  the  variational  equation  (10)  on  each  region 
Furthermore,  the  standard  interface  conditions  will  be  satisfied:  n  x  E  and  n  ■  D  will  be 
continuous  across  interfaces  [2,  pp.  13-16]. 


3.  The  Reflection  Coefficient 

To  facilitate  our  use  of  PMLs  in  the  next  section,  we  compute  the  reflection  coefficient  at 
a  typical  interface  separating  two  diagonally  anisotropic  materials.  For  illustrative  purposes, 
we  do  this  without  loss  of  generality  for  an  interface  consisting  of  the  yz-plane  at  x  —  0. 
For  this  we  generalize  the  calculations  in  [13].  In  the  frequency  domain,  the  time  harmonic 
Maxwell’s  equations  in  a  homogeneous  region  with  no  sources  and  free  charge  density  p  =  0 
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are 

V  •  D  =  0 
V-B  =  0 

(11) 

Vxi?  =  —jtoB 
V  x  H  =  juD  +  Jc. 

The  vectors  in  (11)  are  functions  of  position  f  —  (. x,y,z )  and  frequency  u.  In  the  frequency 
domain  the  constitutive  relations  (3)  for  a  diagonally  anisotropic  material  take  the  form 

D  =  eA  E 

(12)  B  =  HqAH 

Jc  =  a  A  E 


where  e(u)  is  the  complex  permittivity  and 

(13)  A  =  diag(Ax,  Xy,  \z) 

is  a  complex,  frequency  dependent  diagonal  matrix. 

Let 

(14)  ko  =  uy/fjLo(e(u})  + 


Calculations  similar  to  those  in  Section  2  give  the  Helmholtz  equation  for  a  homogenous, 
diagonally  anisotropic  material 

(15)  ((\y\z)  1d‘^c  +  (XxXz)  ldy  +  (XxXy)  ldl)E  +  k^E  =  0. 

Let  A^  =  diag(Aa;,  Xy,  Xz)  be  the  diagonal  matrix  in  the  constitutive  relations  (12)  for 
the  region  x  <  0,  and  let  A+  =  diag^,  jy,  7Z)  be  the  diagonal  matrix  in  the  constitutive 
relations  for  the  region  x  >  0.  Take  the  permittivity  e  and  the  conductivity  a  to  be  the  same 
in  both  regions. 

We  look  for  plane  wave  solutions  of  Maxwell’s  equations  having  the  form 


(16) 


E(t,  r)  =  E0  exp  (— j{k  ■  r  —  cut)) 
H(t,  r)  =  H0  exp  (— j{k  ■  r  —  ut)) 


where  Eq,  H q,  and  k  are  vectors  with  complex  components,  and  j  is  the  imaginary  unit. 
Requiring  that  the  plane  wave  expression  for  E  in  (16)  solve  the  Helmholtz  equation  (15) 
gives  a  dispersion  relation  for  k  =  (. kx ,  ky ,  kz)  in  each  region.  For  x  <  0 


(Aj,A2)  Xk2x  +  (A^A^)  1k2  +  (XxXy)  lk2z  —  k$. 
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For  x  >  0  a  similar  relation  holds  with  7  replacing  A.  The  dispersion  relation  is  the  equation 
of  an  ellipsoid  and  has  solutions 


kx  =  h) \J \y\z  cos  9  cos  0 
ky  =  k0  a/ XxXz  cos  9  sin  0 
kz  =  kQy/XxXy  sin  9 


where  9  and  0  are  complex.  In  our  simulations  E(t,f )  will  be  independent  of  y,  so  we  take 
0  =  0  and  ky  =  0.  Let  x,  y,  and  5  denote  the  standard  basis  vectors.  Then  for  x  <  0  the 
incident  and  reflected  wave  vectors  are 


h  =  k0(^\yXz  cos (9i)x  +  yfXxXy  sin(6>05) 

K  =  kp{—\/xyxz  cos {9r)x  +  \JXxXy  sin(dr)i) 


and  for  x  >  0  the  transmitted  wave  vector  is 


(17) 


kt  =  ^o(y7?/7*  cos (9t)x  +  y/^y  sin(6>t)z) . 


Let  TM.y  denote  a  plane  wave  where  H  has  only  a  ^/-component,  and  let  TEy  denote  a 
plane  wave  where  H  has  only  x-  and  ^-components.  Any  plane  wave  can  be  written  as 
a  superposition  of  TE^  and  TMy  waves.  Let  7 Z  and  T  be  the  reflection  and  transmission 
coefficients  at  the  interface  for  a  TMy  wave.  Then  using  a  complex  parameter  A  we  can 
express  the  incident,  reflected,  and  transmitted  magnetic  fields  as 


Hi  =  y 


h 


0 


Hr  =  y 


Ht=y 


LO/Iq 

ko 

cn/io 

k0 


UflQ 


A  exp (~jk0 (\/XyXz  cos (9i)x  +  \fXxXy  sin(9i)z))  exp(jut) 

7ZA  exp(-jk0(-\/XyXz  cos(9r)x  +  \J XxXy  sin(dr)^))  exp(ju;f) 
T A  exp (~jk0( y/'fy'fz  cos (9t)x  +  y^TT^sin (9t)z))  exp(jujt) 
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Calculating  Vxil  using  Maxwell’s  curl  equations  (11)  and  the  constitutive  relations  (12), 
we  obtain  the  electric  fields 


E,  = 


'  y  sin (9i)x  ~  ]Jy  cos (&i)z 


A  exp  (—jko(y/\y\z  cos  (6i)x  +  \JXxXy  sin  (0i)z))  exp  (jut) 


Er  = 


'  — ^  sin(dr)a:  +  \!  cos(#r)£ 


E A  exp ( —j k0(— \fXyXz  cos  [6r)x  +  yj  XxXy  sin(0r)z))  exp  {jut) 


Et  = 


—  sin(6U£  —  .  —  cos  (9t)z 

lx  V  X* 


TAexp(-jk0(yy/'yy'yzcos(9t)x  +  y/ixiysm(9t)z))  exp  (jut). 

The  tangential  components  of  the  total  electric  and  magnetic  fields  are  continuous  at  the 
interface  x  =  0.  Enforcing  this  continuity,  for  example,  at  the  point  (x,  z)  =  (0,  0)  we  have 

i  +  n  =  r 


(18) 


cos  6i  +  E\  1^-  cos  0r  =  —Tx  j  —  cos  9t. 
A*  V  A, 


Solving  for  E  gives  the  equation  for  the  reflection  coefficient  for  a  TMy  wave: 


(19) 


E  = 


1  Y1  COS  0i-  cos  0t 
1  Y1  cos  0r  +  .  cos  6t 


Taking  the  ratios  of  //-coordinates 


y-  Hr  y  ■  Ht 

- —  and  - — 

V  ■  Hi  y  ■  Hi 

to  be  constant  on  the  interface  x  =  0  gives  Snell’s  law: 

sin  0i  =  sin  9r 


(20) 


V  V.-A,,  -.ill  <>:  =  *JYxXv  sin  8t- 


A  similar  calculation  shows  that  (19)  and  (20)  are  also  the  reflection  coefficient  and  Snell’s 
law  for  TE y  waves  at  the  x  =  0  interface. 


Proposition  3.1.  Assume  that  T  ^  0  and  that 


(21) 


^y/^z  r‘iy/r‘iz 
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Then  equations  (18)  and  (20)  imply  that  1Z  —  0  and  there  exist  integers  n  and  m  such  that 

9r  =  0i  +  2nn 
9t  —  9i  +  2rmr. 

Proof.  By  Snell’s  law  (20)  and  the  Pythagorean  identity  we  must  have  cos  9i  =  ±  cos  9r 
and  cos  9{  =  ±cos  9t-  Since  T^O  the  linear  equations  (18)  will  have  a  finite  solution  if  and 
only  if  cos  9i  =  cos  9r  =  cos  9t.  Consequently  7Z  =  0  by  (19)  and  sin(#j  —  9r)  =  sin (9i  —  9t)  =  0 
by  standard  trigonometric  identities.  The  complex  sine  function  has  only  real  roots,  so  9%  is 
related  to  9r  and  9t  by  (22). 

Thus,  if  A+  and  A-  satisfy  (21),  the  planar  interface  normal  to  the  x-axis  will  be  reflec¬ 
tionless,  and  the  angles  of  incidence,  reflection,  and  transmission  will  be  the  same  for  all 
plane  waves.  Similar  results  can  be  obtained  for  planar  interfaces  normal  to  the  y-axis  or 
z-axis. 


4.  Problem  Formulation 

In  this  section  we  give  a  problem  geometry  that  allows  for  arbitrarily  many  PMLs,  and  for 
each  PML  we  define  the  anisotropy  matrix  T  in  the  constitutive  relations  (3).  We  show  that 
the  PML  interfaces  are  reflectionless  and  that  energy  is  attenuated  in  each  PML.  The  vari¬ 
ational  equation  (10)  on  each  homogeneous  region  is  used  to  formulate  a  global  variational 
equation  for  E  with  the  PML  and  polarization  terms  satisfying  auxiliary  differential  equa¬ 
tions.  Finally,  an  absorbing  boundary  condition  is  placed  on  three  sides  of  our  computation 
domain  to  further  damp  the  reflected  energy. 

Problem  Geometry.  We  assume  that  E,  P ,  and  Js  in  (10)  are  independent  of  y  so  our 
computational  domain  is  contained  in  the  xz- plane  as  depicted  in  Figure  2.  Let  X  and  Z  be 
closed  intervals  in  the  respective  x-  and  z-axes,  and  let  T>  =  X  x  Z  denote  the  computational 
domain. 

We  partition  the  intervals  X  into  disjoint  closed  intervals  AT,  A0,  and  X+  such  that 

max  AT  =  minX0 
max  X0  =  ruin  A+, 

and  partition  the  interval  Z  into  disjoint  closed  intervals  Z_,  Zv,  and  Z j  such  that 

max  2T  =  min  Zv 
max  Zv  =  min  Z(j. 

The  regions  Zv  and  Zci  are  the  vacuum  and  dielectric  regions  respectively.  Let 

Zq  =  Zv  U  Zcj  . 

The  region  X0  x  Z0  will  be  our  domain  of  interest.  The  buffer  region  T>\(Xq  x  Z0)  outside  our 
domain  of  interest  will  contain  the  PMLs.  This  problem  geometry  is  illustrated  in  Figure  2. 
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x-axis 


Antenna 


FIGURE  2.  2-D  problem  geometry 

Perfectly  Matched  Layers.  Let  (3x[x)  and  /3z(z)  be  piecewise  constant,  convex  functions 
such  that  fix  —  0  in  X0  and  (3Z  =  0  in  Z0.  The  functions  fjx  and  (5Z  control  the  damping  in 
the  PML  regions.  The  convexity  requirement  implies  that  fix  >  0,  fiz  >  0  and  that  damping 
will  be  greater  closer  to  the  boundary  of  T>.  Define 

ux(t,x)  =  5(t)  +Px(x)H(t) 
uz{t,z)  =  S(t)  +  (3z{z)H{t) 

where  7d  is  the  Heaviside  function.  Let  Sx  and  Sz  be  convolution  operators  on  C[0,  oo)  with 
kernels  ux  and  uz  respectively.  Then  for  fixed  x  and  z,  Sx  and  Sz  are  bijections  on  C[0,  cxd) 
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since  for  any  real  fi 

(. 8{t )  + 1 3H(t ))  *  (S(t)  -  finfifie^1)  =  S(t). 

Thus  Sx,  Sz  G  B+  with  B+  defined  in  (2). 

On  X  x  Z  let  the  operator  T  in  (3)  be  given  by 

(S?sg  0  0  \ 

(24)  T  —  [  0  SXSZ  0  . 

V  0  0  SXS~  7 

Although  T  is  defined  for  a  two-dimensional  problem  there  is  a  straightforward  generalization 
to  three  dimensions  [9].  Our  domain  of  interest  is  isotropic  since  Sx  and  Sz  are  the  identity 
operator  in  the  region  X0  x  Z0. 

The  operator  T  defined  in  (5)  is 

(25)  f  =  diag(5-2,/,572) 


where  /  is  the  identity  operator. 

We  will  show  that  the  PML  interfaces  are  reflectionless  and  that  the  PMLs  attenuate 
energy.  The  operator  A  in  (12)  corresponding  to  (24)  is 


where 

(26) 


1  T  fix / (jO) 
■s*  =  1  +fiz/(juj). 


Let  x  =  Xq  be  a  point  of  discontinuity  in  fix.  Then  the  interface  x  =  x0  will  be  reflectionless 
by  Proposition  3.1  since  (21)  is  satisfied  when 


A  =  cliag(Ax,  Ay,  A,)  =  lirn  A(x,  z) 

X^XQ 

A+  =  diag(7x,7j/,7J?)  =  lim  A (x,z). 

x— 


Consequently  the  reflection  coefficient  at  the  interface  will  vanish  and  each  plane  wave’s 
angle  of  incidence  equals  its  angle  of  transmission.  Similar  arguments  show  the  interface  at  a 
discontinuity  in  fiz  is  also  reflectionless  with  angle  of  incidence  equaling  angle  of  transmission 
for  incident  plane  waves. 

We  can  show  that  waves  will  be  attenuated  in  each  PML  by  introducing  the  transmitted 
wave  vector  kt  (IT)  into  the  plane  wave  solution  (16).  The  plane  wave’s  exponential  term 
becomes 


(27) 


exp (-jk0(sx(x)  cos (0t)x  +  sz(z) sin(^)e)). 
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In  PMLs  matched  to  free  space  9t  is  real  and  the  value  ko  defined  in  (14)  is  =  ui/c  where 
c  =  1/^eo/io  is  the  wave  speed.  Thus  by  the  definitions  of  sx  and  sz  in  (26)  the  absolute 
value  of  the  exponential  factor  (27)  is 

exp(— (1  /c)(fix(x)  cos (9t)x  +  Pz(z)  sin (9t)z)). 

Thus  in  a  PML  matched  to  free  space  there  is  attenuation  in  the  direction  of  wave  propagation 
cos  9tx  +  sin  9tz,  and  the  rate  of  attenuation  in  the  x  and  z  directions  is  controlled  by  the 
magnitudes  of  j3x  and  /3Z  respectively. 

In  the  dielectric  region  X  x  Zd,  sz  —  1  since  (3Z  —  0  on  Z0,  and  the  product  kocos9t  in 
(27)  is  real-valued  since  planes  of  constant  amplitude  in  the  dielectric  X0  x  Zd  are  parallel 
to  the  xy-plane  [2,  pp.  210-212],  Thus  in  a  PML  matched  to  a  dielectric  the  absolute  value 
of  the  exponential  factor  (27)  is 

exp (—  (Px (x)/u>)k0  cos (9t)x)  exp ((9(/c0  sin  9t))z). 

Thus  there  is  attenuation  in  the  +x  and  —  x  directions  in  the  PMLs  matched  to  a  dielectric, 
the  magnitude  of  /3X  controls  the  rate  of  attenuation,  and  the  attenuation  is  frequency  and 
direction  dependent. 


Auxiliary  Differential  Equations.  We  will  next  reformulate  the  variational  equation  (10) 
so  that  each  PML  term  satisfies  an  auxiliary  differential  equation.  Define  Vx  and  Vz  by 


(28) 


Vx  +  E  —  TXE 
Vz  +  E  =  TZE. 


From  (25)  we  have  Tx  =  Sx2.  So  by  the  definition  of  the  kernel  ux  of  Sx  in  (23)  we  have 

E  =  d*(Tx)~1(Vx  +  E) 

=  d2S2x(Vx  +  E) 

=  d2{ux  *  ux)  *  (Vx  +  E) 

=  (5  +  2pJ{t)+Pl5)*{Vx  +  E) 

=  Vx  +  E  +  2px{Vx  +  E)  +  Pl(Vx  +  E). 

Thus  by  (28)  we  may  replace  the  TXE  term  in  the  variational  equation  (10)  with  Vx  +  E  and 
Vx  satisfies  the  ordinary  differential  equation 

Vx  +  2PX{VX  +  E)  +  P2X(VX  +  E)  =  0. 

Similarly,  since  Tz  =  S~2  we  may  replace  TZE  in  the  variational  equation  (10)  with  Vz  +  E 
where 

Vz  +  2PZ(VZ  +  E)  +  (32z(Vz  +  E)  =  0. 

We  model  the  dielectric  region  X  x  Zd  as  a  Debye  medium  [3,  pp.  11-12],  so  that  the 
polarization  term  P  will  satisfy  the  differential  equation 

(29)  rP  +  P  =  e0edIZdE 
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where  Izd  is  the  indicator  function  on  X  x  Zd ,  r  is  the  relaxation  time  constant,  and  ed  is  a 
constant  measuring  the  difference  between  the  static  relative  permittivity  (when  u  =  0)  and 
the  relative  permittivity  er. 

From  the  preceding  discussion  it  follows  that  the  variational  equation  (10)  in  any  region 
where  f3x  and  /3Z  are  constant  can  be  written  as 

(30)  fioeo(erE,  0)  +  [Iq{P ,  0)  +  /R)( >  0)  +  ( dx(E  +  Vx),  dx(j))  +  ( dz(E  +  Vz),  <920) 

{9X  (E  +  Vx")  j  Ri0) (-£/-[-  Vz) ,  U20) /^0  (4s,  0)  • 

where  the  polarization  and  PML  terms  in  (30)  satisfy 

rP  +  P  =  €0€dIzdE 

(31)  V4  +  2/3,14  +  ftVx  =  -2(3XE  -  ftE 

Vz  +  2/3,14  +  &V*  =  —2{3ZE  -  (3 2zE. 

The  zero  initial  conditions  for  our  problem  imply  that  14  will  be  supported  in  X  x  Z_,VX 
will  be  supported  in  (X_  x  Z)  U  (X+  x  Z),  and  P  will  be  supported  in  X  x  Zd. 

Global  Variational  Equation.  Our  source  term  Js  in  (8)  will  model  an  antenna  in  free 
space,  and  we  assume  the  signal  is  polarized  with  oscillations  in  the  £2-plane  only.  Let 

Js(t,x,z)  =  I{0,tf)(t)I{-x,x)(x)5(z)sm(ut)x 

\6Z)  T  ( .  \  ^ 

=  Js{t,  x,  z)x 

where  I(a,b)  is  the  indicator  function  on  the  interval  (a,  b),  (±x,  0)  G  X0  x  Zv ,  utf  is  an 
integral  multiple  of  27r,  and  Js  is  defined  in  (9). 

There  are  interfaces  between  PMLs  at  jumps  in  the  PML  parameters  (3x{x)  and  /3z(z). 
By  construction  these  PML  interfaces  are  reflectionless.  However,  when  the  problem  is 
discretized  numerical  reflections  occur  at  PML  interfaces.  These  numerical  reflections  can 
be  made  arbitrarily  small  by  making  the  jumps  in  /3X  and  /3Z  sufficiently  small. 

We  note  that  during  a  finite  time  interval  [0, 4nai],  the  sum  of  the  boundary  integrals 
in  (30),  {dx(E  +  Vx),  nx(j))dn  or  {dz(E  +  Vz),  nz4>)dni  that  arise  from  both  sides  of  a  PML 
interface  can  be  made  arbitrarily  small  by  taking  the  jumps  in  (3X  and  f3z  to  be  sufficiently 
small.  Comparison  calculations  we  have  performed  with  and  without  the  boundary  integrals 
from  the  PML  interfaces  (i.e.,  numerical  computations  including  the  integral  terms  properly 
with  quadratures  vs.  approximating  the  net  integral  contributions  by  zero)  reveal  that  these 
interface  integral  terms  contribute  very  little  to  understanding  the  error  control  issues  related 
to  reflections  at  the  interfaces  of  the  dielectric  relative  to  those  at  the  boundaries  of  the  finite 
computational  domain.  We  have  therefore,  in  this  paper,  chosen  to  approximate  the  interface 
integrals  by  zero  in  our  calculations,  simplifying  and  concentrating  on  the  reflections  from 
the  computational  boundary  and  the  vacuum-dielectric  interface.  Other  calculations  (not 
reported  here)  verify  that  the  conclusions  of  our  results  with  respect  to  inverse  problem 
methodology  are  not  altered  by  inclusion  of  the  PML  interface  integral  terms  in  our  system. 
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After  approximating  the  sum  of  the  boundary  integrals  from  both  sides  of  a  PML  interface 
by  zero  we  obtain  a  variational  equation  that  will  hold  on  the  entire  domain  V : 

(33)  ho€q (erE,  0)  +  fioiP ,  0)  +  fio(aE ,  0)  +  ( dx(E  +  Vx),  dx<j>)  +  ( dz(E  +  14),  <920) 

(&x(E  +  14),  Tlx(f>') QT>  (9z(E  +  14),  U20)g-p  /i(j ( j.s  1  0)- 

where  the  polarization  and  PML  terms  satisfy  (31). 

Absorbing  Boundary  Condition.  We  next  add  boundary  conditions  to  our  problem 
formulation.  Let  C-x,  C+x,  C_2,  and  C+z  denote  the  four  boundaries  of  V  with  outward 
normals  in  the  —x,  +x,  —z  and  +z  directions  respectively.  Assign  the  boundary  conditions 

(A*  +  dt)E  =  cdxE  on  C_x 

(A*  +  9t)E  =  - cdxE  on  C+x 

(Az  +  dt)E  =  cdzE  on  C-z 

E  =  0  on  C+z 

where  c  =  I/a/zioCo-  The  boundary  condition  on  C+z  models  a  dielectric  with  a  supraconduc- 
tive  backing.  We  show  below  that  the  remaining  three  boundary  conditions  are  absorbing 
boundary  conditions  for  PMLs  matched  to  free  space. 

Take  an  incident  plane  wave  on  the  C+x  boundary  with  the  form 

Ei  =  exp(—  jko(sx(x)  cos {6i)x  +  sz(z )  sin(6)j)^))  exp(jut). 

Then  the  reflected  wave  will  have  the  form 

Er  =  lZex.p(—jk0(—sx(x)  cos(9i)x  +  sz(z )  sin(6lj)^))  exp(jcat). 

The  total  wave  is  E  =  Ei  +  Er. 

For  PMLs  matched  to  free  space,  6i  is  real-valued  and 

kocos(6i)  =  (c o/c)  cos(6)0 

by  (14).  For  PMLs  matched  to  the  dielectric,  we  may  assume  the  plane  wave  traveled  across 
the  vacuum-dielectric  interface.  Applying  Snell’s  law  at  the  vacuum-dielectric  interface 
shows  that  ko  cos (6**)  is  real-valued  [2,  pp.  210-212]  and 

ko  cos(9i)  =  (c o/c)  cos (9i) 

where  9,L  is  the  real-valued  angle  of  incidence  on  the  vacuum-dielectric  interface.  Let  9  =  9t 
for  PMLs  matched  to  free  space  and  9  =  9{  for  PMLs  matched  to  the  dielectric.  Then 
applying  the  boundary  condition  (/3X  +  dt)E  =  —cdxE  on  C+x  leads  to  the  equality 

(Px+ju){  1  +  E)  =  c(jsxuj/c)  cos(0)(l  -  TV) 

=  ju(  1  +  Px/(juj))  cos(6»)(l  -  7 Z) 

=  ((3x+ju)  cos(0)(l  -TZ). 


Solving  for  7 Z  gives 
(35) 
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cos  9—1 
cos  9  +  1 

where  9  is  a  real-valued  angle  of  incidence  for  either  the  boundary  C_x  or  the  vacuum- 
dielectric  interface.  Similarly  one  can  also  show  that  (35)  is  the  reflection  coefficient  for  the 
boundary  conditions  (34)  on  C_x  and  C-z. 

Thus  the  reflection  coefficient  1Z  depends  on  an  angle  of  incidence  and  is  independent 
of  the  frequency.  The  reflection  coefficient  vanishes  for  9  =  0  and  \TZ\  is  an  increasing 
function  of  9  on  (0,  7t/2).  In  our  problem  we  can  generally  expect  9  <  7t/4.  Reflected  waves 
with  sufficiently  large  angles  of  incidence  will  reflect  off  two  boundaries  before  returning  to 
the  domain  of  interest.  Consequently,  it  is  reasonable  to  expect  this  absorbing  boundary 
condition  will  reflect  approximately 

c°s(x/4)  -  _  003 

cos(7t/4)  +  1/ 

of  the  energy.  This  shows  that  the  absorbing  boundary  condition  plays  a  useful  role  in  damp¬ 
ing  outgoing  energy.  Indeed,  unlike  the  situation  suggested  in  [12],  the  absorbing  boundary 
conditions  combined  with  the  PMLs  can  enhance  energy  absorption.  Our  comparison  calcu¬ 
lations  have  verified  this,  at  least  for  the  problems  considered  in  this  paper. 

Since  differentiation  commutes  with  convolution,  it  follows  from  the  definitions  of  Vx  and 
Vz  in  (28)  that  the  same  boundary  conditions  will  hold  with  Vx  and  Vz  replacing  E.  Thus 
the  variational  equation  (33)  becomes 

(36)  Hoeo(erE,  0)  +  Ho{P ,  0)  +  Ho(&B,  0)  +  ( dx(E  +  Vx ),  dx(ft)  +  ( dz{E  +  Vz ),  <920 ) 

+  ((l/c)(Ar  +  dt)(E  +  Vx),4>)c_xuc+X  +  ((Vc)(As  +  dt)(E  +  Vz),4>)c_z  =  —fJ>o  (4,0)- 

5.  Discretization 

To  aid  the  numerical  calculations  we  scale  the  time  variable  by  a  factor  of  c  =  1/^eoJP) 
so  that  E  is  replaced  by  cE,  P  is  replaced  by  cP,  etc.,  and  we  scale  the  polarization  P  by  a 
factor  of  1/co  so  that  P  is  replaced  by  €qP.  The  variational  equation  (36)  becomes 

(37)  ( erE ,  0)  +  (P ,  0)  +  T)o(aE ,  0)  +  (dx(E  +  Vx ),  <9X0)  +  (d Z{E  +  Vz ),  <920) 

+  {{olx  +  dt)(E  +  14),  4>)c-xuc+x  +  ((«z  +  dt)(E  +  14),  0)c_2  =  —Vo  (4, 0) 

where  770  =  \J /i0 /+  is  the  impedance  of  free  space,  ax  =  f3x/c,  and  az  =  f3z/c. 

The  differential  equations  (31)  become 

P  +  +P  =  td,vIzdE 

Vx  +  2otxVx  +  a2xVx  =  —2  (xxE  —  axE 
Vz  +  2azVz  +  oq  14  =  — 2azE  —  o+E 


(38) 
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with  7  =  (cr)  l. 

After  scaling  our  source  term  (32)  in  the  time  variable,  we  obtain 
(39)  Js(t,x,z )  =  I^tf)(t/c)I{-x,x)(x)6(z)  sm((u /c)t) 

where  (u/c)tf  is  an  integral  multiple  of  27 r. 

We  will  use  a  Galerkin  finite  element  approximation  to  discretize  the  problem  in  the 
spatial  variable  to  obtain  piecewise  linear  approximates  for  E(t,  •)  and  P(t,  •).  Choose  strictly 
increasing  sequences  of  real  numbers 


ii,  =  ujy  c  x 
n,  =  c  z 

such  that  n.r  contains  the  end  points  of  X  and  such  that  hh  contains  the  end  points  of  Z 
and  Zd.  Our  grid  will  be  the  Cartesian  product 

(40)  Q  =  {(x,z)  :  x  G  n^,  z  G  IIZ}, 

and  the  vacuum-dielectric  boundary  will  coincide  with  a  line  of  grid  points.  We  then  con¬ 
struct  a  triangulation  T  of  the  computational  region  T>  such  that  for  each  triangle  r  G  T, 
the  set  of  vertices  v{r)  satisfies 


V(T)  =  {{Xi,Zj),  (Xi,Zj+ 1),  (Xj+i,  Zj+i)} 


or 


v{t)  =  {( Xi,Zj ),  (Xi+i ,Zj),  (xi+1,zj+1)} 


for  some  1  <  i  <  Nx  and  1  <  j  <  Nz. 

Since  E  =  0  on  C+z,  let  {qj}f=1  be  an  enumeration  of  the  N  =  NXN~  —  (Nx  —  2)  grid  points 
lying  in  D  that  are  not  in  the  interior  C+z  and  let  {(pi(x,  z)}^=1  be  the  standard  piecewise 
linear  spline  functions  that  are  linear  on  each  triangle  r  G  T  and  satisfy 


MOj)  = 


1  i  —  j 

0  i^j. 
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Then  we  can  approximate  E,  P,  Vx ,  and  Vz  with 


(41) 


N 

E(t,x,z)  fa  ej(t)(f>j(x,  z) 

3  =  1 
N 

P(t,x,z )  fa 

3=1 

N 

Vx(t,x,z)  ^ 

3=1 

N 

Vz(t,x,z)  fa  ^(v2)j(t)^(x,z). 

3=1 


The  basis  functions  4>j  and  the  coefficient  functions  ej,  pj ,  (vx)j,  and  (vz)j  depend  on  the 
chosen  grid  Q.  Let  e,  p,  vx,  and  vz  be  TV-dimensional  vectors  formed  from  these  coefficient 
functions  (e.g.,  e  =  (ei,  . . . ,  ejv)T). 

If  <f>  =  cj4>j i  then  we  will  approximate  ax(x)(f>(x,  z )  and  az(z)(f)(x,  z)  with 


(42) 


N 

olA  fa  ax{x{qj))cj4>3 

3=1 

N 

aA  fa  ^^az(z(qj))cj(j)j. 

3=1 


where  x(qj)  and  z{q3)  denote  the  x  and  z  coordinates  of  qj. 

Making  the  substitution  <f>  =  fa  and  using  the  approximations  (41)  and  (42)  in  the 
variational  equation  (37)  and  the  differential  equations  (38),  we  obtain  a  Galerkin  finite¬ 
dimensional  system  of  equations 


(43)  (M  +  MzAoo  ~~  l))e  +  Mzdp  +  (Lx  +  AxBx)vx  +  (. Lz  +  AzBz)vz 

+  (Lx  +  Lz  +  AXBX  +  AzBz)e  +  Bxvx  +  Bzvz  +  ( Bx  +  Bz  +  rjo  crMzd)e  =  J 

with  the  polarization  and  PML  terms  satisfying 

P  +  IP  =  tdlhe 

vx  +  2 Axvx  +  A2xvx  =  -2 Axe  -  A2xe 
vz  +  2 Azvz  +  A2zvz  =  -2 Aze  -  A2ze 


(44) 
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where  is  the  relative  permittivity  in  the  dielectric  and 


(45) 


Mij 

(Mzd)ij 

(Lx)i,j 

( BX)i,j 
( BZ)i,j 


(A 


x)i,3 


Ji 


{(ftj,  4>i) 

( izd4>jAi ) 

(9x4*3 1  9x4*i) 

(dz(pj i  9z<t>i) 

(4>j,4>i)c-xuc+x 

(<t>j,<t>i)c-z 

f  ax(x(qj))  i  —  j 

\  0  i^j 

f  az(z(qj))  i  =  j 

\  0  i  ^  j 

f  1  i  =  j  and  IZd4>i  ^  0 
1  0  otherwise 

Vo(jsi  4*i)  ■ 


We  replace  the  MZdp  term  in  (43)  by  using  the  differential  equation  and  its  derivative  for 
P  in  (44)  to  obtain 

MzJj  =  rfMZdp  +  ed7 MZde  -  ecrf2MZde. 

The  hnal  form  of  our  variational  equation  is 


(46)  (M  +  MZd(e  oo  —  l))e  +  +  (Lx  +  AxBx)vx  +  (Lz  +  AzBz)vz 

+  ( Lx  +  Lz  +  AXBX  +  AZBZ  —  ed'y2  MZd)e  +  Bxvx  +  Bzvz 

+  (Bx  +  Bz  +  ( rjacr  —  e,fi)MZdl)e  =  J 

with  the  polarization  and  PML  terms  satisfying  (44). 

The  resulting  system  of  equations  can  be  written  as  a  first  order  system  in  the  composite 
variable 

i  =  (p,vx,vz,e,vx,vz,e)T, 

and  we  approximate  a  solution  to  this  system  using  the  Crank-Nicholson  scheme  [3,  pp. 
58-60], 


6.  Numerical  Results 

In  this  section  we  present  numerical  results  which  show  that  PMLs  attenuate  energy  and 
that  numerical  reflections  from  the  PML  interfaces  can  be  controlled. 

For  these  simulations  the  conductivity  is  cr  =  1.0  x  10~2,  and  the  Debye  parameters  in  the 
differential  equation  (29)  are  r  =  1.0  x  lO^11  and  =  30.0.  The  relative  permittivity  in  the 
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dielectric  is  Coo  =  5.0.  The  antenna  is  centered  at  the  origin,  and  the  interrogating  signal 
(39)  will  have  3  cycles.  The  source  term  parameters  in  (39)  are 

x  =  0.03 

u  —  2n  x  1.8  x  109  rad/sec 
tf  =  67Tc/cu  (~  0.5  sec). 

The  time  step  for  the  finite  differencing  is  1.0  x  10“4  units  in  the  scaled  time  variable. 

We  define  the  following  sequence  to  discretize  our  PML  regions:  £o  =  0.05  and 


0  =  0-1  +  (0-0025)  +  ^(0.0125) 

j  —  1,  2, ...  60 

For  the  PML  parameters  let  cto  =  0.0  and 

(47)  a,  =  +  (0.01)  +  ^(0.11) 

J  =  1,2,...  60. 

Let 

x_  =  k60,-O] 

Z-  —  [— Oo,  —O] 

X0  =  [-0.05,  0.05] 

Zv  =  [-0.05,0.05] 

x+  =  [£o,6o] 

Zd  =  [0.05,0.35], 

Our  grid  Q  will  be  the  Cartesian  product 

q  —  {(x,  z) :  x  g  n  u  n„  u  n+, 

zeLu  n„  u  nfi} 

where 


n„  =  {(0.0025)  j  :  -20  <j<  20} 

Ud  =  {(0.0025)  j  :  20  <j<  140} 

n+  =  (0  :  0  <  j  <  60} 
n_  =  (-0  :  o  <j<  60}. 

Let  0^(0  be  any  piecewise  constant,  convex  function  satisfying  ax  —  0  on  X0  and  ax(Cj)  = 
ax(—£,j)  =  otj.  Let  az  be  any  piecewise  constant,  convex  function  satisfying  az  =  0  on 
Z0  =  Zv  IS  Zd  and  (-£,■)  =  «r 

There  will  be  a  loss  in  accuracy  clue  to  the  coarser  grid  as  we  approach  the  boundary  of 
D  =  X  x  Z ,  but  the  coarser  grid  occurs  outside  of  our  domain  of  interest  and  substantially 
reduces  the  number  of  grid  nodes  needed  for  our  simulation. 

The  interfaces  between  PMLs  are  reflectionless;  however,  when  the  problem  is  discretized 
numerical  reflections  occur  when  there  are  jumps  in  the  PML  parameter.  These  numerical 
reflections  at  a  PML  interface  may  be  made  arbitrarily  small  by  making  the  jump  in  the 
PML  parameters  oq  sufficiently  small. 

We  have  also  computed  the  solution  when  the  grid  Q  is  preserved  and  the  PMLs  are 
replaced  with  free  space.  The  absolute  difference  between  the  free  space  solution  and  the 
PML  solution  in  the  domain  of  interest  A"o  x  Z$  is  less  than  0.6.  By  contrast,  the  amplitude  of 
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the  signal  after  reflecting  off  the  vacuum-dielectric  interface  is  approximately  15  as  we  will  see 
in  Figure  3.  This  implies  that  the  error  in  the  domain  of  interest  due  to  numerical  reflections 
from  the  PMLs  is  less  than  0.6/15  =  4%  of  the  amplitude  of  the  reflected  interrogating 
signal. 

In  Figures  3,  4,  and  5  we  give  snapshots  E(t,0,z)  of  the  electric  held  along  the  2-axis  at 
times  t  =  0.6001,  1.6001,  and  3.6000,  respectively. 

The  source  current  Js  propagates  three  cycles  during  the  time  interval  0  <  t  <  tf  ~  0.5, 
and  energy  travels  at  unit  speed  in  the  vacuum  region  —0.05  <  2  <  0.05  and  the  PML  region 
2  <  —0.5.  Consequently  at  time  t  =  0.6001  (Figure  3)  part  of  the  original  signal  (the  left 
moving  part:  the  antenna  emits  both  left  and  right  moving  interrogating  signals  into  the 
vacuum)  will  be  in  the  2-interval  (—0.6,  —0.1),  the  signal  reflected  from  the  vacuum-dielectric 
interface  at  2  =  0.05  will  be  in  the  2-interval  (—0.5,  0),  and  the  transmitted  signal  will  be  in 
the  dielectric  region  0.05  <  2  <  .35.  The  trailing  wave  crest  of  the  reflected  signal  is  visible 
at  2  =  0  and  has  an  amplitude  of  approximately  15.  The  PML  induced  damping  of  the 
signal  is  evident  in  the  PML  region  2  <  —0.05. 

At  time  t  =  1.6001  (Figure  4)  the  Brillouin  precursors  [1,  3,  8]  can  be  seen  in  the  Debye 
medium  0.05  <  2  <  0.35.  The  signal  in  the  Debye  medium  is  still  moving  in  the  positive  2 
direction.  The  original  left-moving  signal  and  the  signal  reflected  from  the  vacuum-dielectric 
interface  have  had  sufficient  time  to  hit  the  computational  boundary  and  return  to  the 
domain  of  interest  A"o  x  Zq.  All  the  energy  evident  in  the  vacuum  interval  —0.05  <  2  < 
0.05  has  returned  after  reflecting  off  a  PML  interface  or  the  computational  boundary.  The 
amplitude  of  this  reflected  energy  is  less  than  0.4  while  the  amplitude  of  the  leading  wave 
crest  in  the  dielectric  is  approximately  1.6. 

At  time  t  =  3.6  (Figure  5)  the  transmitted  signal  has  reflected  off  the  dielectric’s  supra- 
conductive  backing  and  has  emerged  from  the  dielectric.  The  leading  wave  crest  of  this 
secondary  reflection  is  visible  at  2  =  0.  The  energy  that  has  returned  to  the  domain  of 
interest  Xq  x  Zq  after  reflecting  off  a  PML  interface  or  the  computational  boundary  is  seen 
as  “saw  teeth”  superposed  on  the  leading  wave  crest  of  the  secondary  reflection  off  the  dielec¬ 
tric’s  supraconductive  backing.  This  figure  shows  that  the  depth  of  the  dielectric  material  is 
recoverable  in  the  inverse  problem  (see  [3]  for  details  on  how  this  is  done). 

The  amount  of  numerical  reflection  from  the  PML  interfaces  can  be  controlled  by  making 
the  jumps  in  the  PML  parameters  sufficiently  small.  For  instance,  if  the  PML  parameters 
Oj  in  (47)  are  redefined  to  be 

oij  =  otj_ i  +  (0.001)  +  ^——^(0.011)  j  =  1,  2, . . .  60. 

59 

then  a  simulation  shows  that  the  error  in  the  domain  of  interest  due  to  numerical  reflections 
from  the  PMLs  is  less  than  0.08.  (As  mentioned  above,  the  PML  parameters  in  (47)  give  a 
maximum  error  of  0.6.)  However,  the  damping  in  the  PML  region  is  much  less  with  these 
smaller  parameters.  If  smaller  jumps  in  the  PML  parameter  are  used  then  more  grid  points 
are  needed  in  the  regions  AT,  X+,  and  to  obtain  the  damping  achieved  using  larger 
jumps. 
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Figure  5.  Graph  of  E( 3.6,  0,  z)  for  z  G  Z 


7.  Concluding  Remarks 

Numerical  simulations  demonstrate  that  reflections  from  interfaces  at  PML  boundaries 
can  be  controlled  by  reducing  the  number  or  magnitude  of  the  jumps  in  the  PML  parameter. 
Furthermore,  our  simulations  have  demonstrated  that  attenuation  occurs  in  the  PML  regions, 
so  expanding  the  size  of  PML  regions  results  in  greater  energy  attenuation.  Thus  a  strategy 
of  reducing  the  size  and  number  of  jumps  in  the  PML  parameter  while  expanding  the  size 
of  the  PML  regions  will  allow  one  to  simultaneously  control  both  the  reflection  from  PML 
interfaces  and  the  return  of  energy  from  the  computational  boundary. 

The  successful  use  of  PMLs  in  these  simulations  enables  us  to  compute  the  electromagnetic 
field  and  solve  a  2-D  forward  problem  while  controlling  reflections  from  finite  computational 
domains.  We  fully  anticipate  that  our  solution  methodology  for  the  forward  problem  will 
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permit  us  to  determine  dielectric  depth  and  parameters  in  an  inverse  problem  thereby  ex¬ 
tending  the  inverse  problem  methodology  for  the  one-dimensional  case  given  in  [3].  We  are 
currently  combining  the  ideas  presented  here  with  a  reduced  order  computational  method¬ 
ology  (Proper  Orthogonal  Decomposition  based  techniques,  see  [4]  for  details)  in  pursuing 
multidimensional  domain  inverse  problems  that  generalize  to  higher  dimensions  the  ideas  of 
[3], 


8.  Acknowledgments 

The  authors  are  grateful  to  Richard  Albanese  for  numerous  discussions  and  frequent  en¬ 
couragement.  This  research  was  partially  supported  by  the  U.  S.  Air  Force  Office  of  Scientific 
Research  under  grant  AFOSR  F49620-01-1-0026. 

References 

[1]  Richard  Albanese,  John  Penn,  and  Richard  Medina.  Short-rise-time  microwave  pulse  propagation 
through  dispersive  biological  media.  J.  Optical  Society  of  America  A,  6(9):1441-1446,  1989. 

[2]  Constantine  A.  Balanis.  Advanced  Engineering  Electromagnetics.  John  Wiley  &  Sons,  New  York,  1989. 

[3]  H.  T.  Banks,  M.  W.  Buksas,  and  T.  Lin.  Electromagnetic  Material  Interrogation  Using  Conductive 
Interfaces  and  Acoustic  Wavefronts.  SIAM,  Philadelphia,  2000. 

[4]  H.  T.  Banks  and  G.  M.  Kepler.  Reduced  order  computational  methods  for  electromagnetic  material 
interrogation  using  pulsed  signals  and  conductive  reflecting  interfaces.  J.  Inverse  and  Ill-posed  Problems , 
11(4),  2003. 

[5]  Jean-Pierre  Berenger.  A  perfectly  matched  layer  for  the  absorption  of  electromagnetic  waves.  J.  Comput. 
Physics ,  114(2)  :185— 200,  1994. 

[6]  J.  G.  Blaschak  and  J.  Franzen.  Precursor  propagation  in  dispersive  media  from  short-rise-tinre  pulses 
at  oblique  incidence.  J.  Optical  Society  of  America  A ,  12:1501-1512,  1995. 

[7]  V.  A.  Bokil  and  M.  W.  Buksas.  A  2D  mixed  finite  element  formulation  of  the  uniaxial  perfectly  matched 
layer.  J.  Comput.  Physics,  submitted. 

[8]  L.  Brillouin.  Wave  Propagation  and  Group  Velocity.  Academic  Press,  New  York,  1960. 

[9]  M.  W.  Buksas.  Implementing  the  perfectly  matched  layer  absorbing  boundary  condition  with  mimetic 
differenceing  schemes.  In  F.  L.  Teixeira,  editor,  Geometric  Methods  for  Computational  Electromagnetics , 
number  32  in  Progress  in  Electromagnetic  Research  Series,  pages  383-411.  EMW  Publishing,  Cambridge, 
Mass.,  2001. 

[10]  Bjorn  Engquist  and  Andrew  Majda.  Absorbing  boundary  conditions  for  the  numerical  simulation  of 
waves.  Mathematics  of  Computation,  31(139):629-651,  1977. 

[11]  F.  G.  Friedlander.  Introduction  to  the  Theory  of  Distributions.  Cambridge  University  Press,  Cambridge, 
1982. 

[12]  P.  G.  Petropoulos.  On  the  termination  of  the  perfectly  matched  layer  with  local  absorbing  boundary 
conditions.  J.  Comput.  Physics,  143:665-673,  1998. 

[13]  Zachary  S.  Sacks,  David  M.  Kingsland,  Robert  Lee,  and  Jim-Fee  Lee.  A  perfectly  matched  anisotropic 
absorber  for  use  as  an  absorbing  boundary  condition.  IEEE  Trans.  Antennas  Propagat .,  43(12):1460- 
1463,  December  1995. 

[14]  E.  Turkcl  and  A.  Yefet.  Absorbing  PML  boundary  layers  for  wave-like  equations.  Appl.  Num.  Math., 
27:533-557,  1998. 

[15]  R.  W.  Ziolkowski.  Time-derivative  Lorentz  material  model-based  absorbing  boundary  condition.  IEEE 
Trans.  Antennas  Propagation,  45:1530-1535,  1997. 


TIME  DOMAIN  ELECTROMAGNETIC  SCATTERING 


27 


Current  address,  Brian  L.  Browning:  GlaxoSmithKline,  Five  Moore  Drive,  Research  Triangle  Park,  NC 
27709,  USA 


